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TECHNICAL MEMORANDUM X- 53365 


BOUNDARY VALUE PROBLEMS ASSOCIATED WITH OPTIMIZATION THEORY 

SUMMARY 


Optimization theory is applied to the physics of trajectory problems to 
yield a boundary value formulation. Three practical numerical methods for 
obtaining solutions to boundary value problems are discussed, and a detailed 
explanation of tlie application of the most efficient method to a sample problem 
is included. Numerical methods are needed because the application of optimi- 
zation theory to trajectory problems generally results in nonlinear differential 
equations connecting the boundary conditions that are not all specified at the 
same value of the independent variable. The methods discussed are adaptable 
to the solution of boundary value problems other than the ones associated with 
trajectory problems or optimization theory. 


INTRODUCTION 


The calculation of trajectories or flight paths for rocket-powered vehicles 
and projectiles has long been of interest to mathematicians, physicists, and 
engineers. With the advent of modern steerable rocket-powered vehicles the 
problem has become especially interesting because of the multiplicity of paths 
that can occur. Obviously, for a particular vehicle, there should be some path 
that allows the attainment of a desired destination most economically. More 
generally, the vehicle must be steered so that it achieves a certain mission in 
an optimum fashion, where optimum usually means most economically although 
other considerations may be added to the definition. The determination of an 
optimum path can be accomplished by the application of optimization theory 
to the physics of the problem. The next section of this report will be concerned 
with an explanation of this procedure with the aid of a simple example. 


GENERAL DISCUSSION 


Formulation of the Boundary Value Problem 


From elementary physics it is known that a force acting on a point mass 
produces an acceleration. Locate the force and point mass in a two-dimensional 
Cartesian coordinate system with its origin at the center of a gravitational field 
and the most simple mathematical approximation to the motion of a rocket- 
powered vehicle is obtained. To be more specific the following diagram is used. 


y y 



The diagram aids in the expression of the differential equations of motion, 
x = (F/m) sin x - g 

X 

y = (F/m) cos x ~ g 

The equations are the sums of the accelerations in the x and y directions; no 
aerodynamic accelerations are considered. F is assumed to remain constant 
for the duration of the powered flight, and the mass of the rocket is assumed to 
diminish at a constant rate (i. e.,fuel is burned at a constant rate) . Then, 


2 


m(t) = m 0 - m(t - t 0 ) 



i 



where m is the constant rate of change of mass with respect to time. Also, 
the accelerations caused by gravity can be expressed as functions of the position 
of the vehicle because gravity fields are usually assumed to be potential fields. 
The potential per unit mass outside a large spherical mass is given by the ex- 
pression 
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y and ^ is a constant characteristic of a particular gravity field. 


r~r 

where R =\ 

r Pl^ v* 4 - V* V\ t r f 4-1 rt 1 f 1 ^>1 /4 
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The preceding relationships define all of the expressions in the equations 
for x and y except the angle x* Chi (x) is called the steering angle or the control 
angle, and it is to be determined at every instant of the powered flight so that 
the thrust (F) will be pointed in an optimum direction to minimize some specified 
performance parameter. If a steering angle program x(t) and starting condi- 
tions (xd, y 0 , x 0 , y„, t 0 ) are given, then the equations of motion (x and y ) can 
be integrated numerically to describe the flight path; and x(t) , y(t) , x(t) ,y(t)are 
determined and tabulated at intervals of time. The symbols x and y are the first 
integrals of the equations of motion. 

Because many different steering angle programs could be specified for 
this example, the problem is to find a particular program that minimizes the 
fuel consumed to move the vehicle from the starting conditions to the desired 
end conditions. For example, if the mission is to achieve a circular orbit at a 
specified altitude (R^) , the following conditions must be satisfied at the termi- 
nation of the thrusting period for the rocket-powered vehicle: 




R f x 2 + y 2 


# 


v f =*I^ + y 2 =J 

JL. 

R f 



0 = xx + yy 

t 


or 

Gj = x 2 + y 2 - R f 2 

= 0 

; 


G 2 = x 2 + y 2 - v f 2 

= 0 



G 3 = xx + yy 

= 0 . 
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The second expression (v = ) is determined by equating the force caused 

I A J Rj 

by gravity to the centrifugal force, due to the rocket's velocity. The third ex- 
pression states that the velocity vector must be perpendicular to the radius 
vector (i. e. , their dot product is zero) . 

A more abstract notation will now be used for convenience in explaining 
the conditions that must be satisfied by the optimum steering angle program. 

The differential equations of motion are written in the following first-order form: 

x. = f.(x , u, t) 

where 

i = 1,2, • • • ,n 

x= Xj, x 2 , • » • ,x (state variables) 
u= u 1( u 2 , • • • , u (control variables) . 

For the example problem the substitutions required are: 

Xj = x 

x 2 =y 

x 3 = X 

x 4 = y 

x 5 = m 
and Uj = x 

Then x t = x s 
x 2 =x 4 



x 5 = - m 


The quantity to be optimized (maximized or minimized) is written as follows: 
n 

s = Z c i x i ( V 

i=l 

where the c^'s are arbitrary constants and t^ is the time at which the thrust is 

terminated and the mission conditions are to be satisfied. 
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For the example problem, choose 


c ! = c 2 = c 3 = c 4 = 0, and c 5 = 1 . 

Then S = m(t^) and a maximization of S will cause the fuel used to be a minimum. 
Also, the starting conditions and mission conditions are written in a functional 
form: 

F (x°, t°) = 0 , 
a 

where x° = ( Xi°, x 2 °, • • ♦ , x^ 0 ) 
a = 1,2, • • • ,k ^ n + 1 - 

WV 0 - 

_ f f f 

where x^ = ( x 4 , x 2 , • • ♦ ,x r ) 

/3 = 1,2, • • • ,1 < n+ 1 

For the example system of differential equations, the functional form of 
the boundary conditions can be written explicitly. 

F i = x i (V - x 0 = 0 

F 2 = *2 (to) -y 0 = 0 

F 3 = X 3 (to) - X 0 = 0 
f 4 = *4 (to) - y 0 = 0 
F s = xg (to) - m 0 = 0 
F e = t (W -t 0 = 0 

Gj= (x 2 4-y 2 ) -R 2 ^ = 0 
G 2 = (x 2 + y 2 ) - v2(t f ) - 0 
G 3 = xx + yy = 0 . 
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Now the flight path optimization problem can be stated completely in the 
more abstract notation. The statement of the problem is as follows: Determine 
a TT(t) that satisfies the system of differential equations 


x. = f. (x, u, t) 
11 


with boundary conditions 


i = 1,2, • • • ,n 


x =x t , x 2 ,» •• ,x n 
a = u t, u 2 , • •• ,« r 


( x° , t°) = 0 a = 1 , 2 , ♦ • • , 1 < n + 1 

G /3 (*f * V = 0 /3 = 1 , 2, • • • , m < n + 1 

and that maximizes or minimizes the quantity 
n 

S= Z c i x i<V • 

i=l 


To determine a u(t) that satisfies the problem statement, an adjoint 
system of differential equations must be defined 


n 


w = - Yj x. 
1 i=i J 


3f. (x, u, t) 

J — ^ 

3 x. 

l 


i = 1 , 2, • • • ,n 
j = 1, 2, • • • ,n 


Now the total system of differential equations (x. and X .) can be integra- 
ed numerically to yield x^(t) and ?y(t) if the initial conditions 1 (x.°, X^°, and t°) 

are known and a steering function u^(t) is specified in some manner. The 
optimum steering functions ^(t) can be specified by defining 
n 

H = 2 A. f.(x, u,t) . 

j=l 3 3 


Then another necessary condition to be satisfied for a minimization 
(maximization) of the quantity S is that the function H be a maximum (minimum) 
with respect to u at every to — t ^ t^. Actually, the maximization or minimiza- 

3 H 

tion of H usually allows the relations - — =0 (k = 1,2, . . . ,r) to be solved to 

3 u, 
k 

give a steering function u. that depends on the x.'s and A.’s at every t 0 — t ^ t . 

K 11 1 
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Also, the second order terms in the series expansion for H in terms of about 

the that makes the first partials zero must be examined to assure a maximi- 
zation or minimization of H with respect to u^. This examination results in a 

condition which requires that the eigenvalues or the principle minors of the 
matrix of second partial derivatives of H with respect to u^ be all positive for 

H a minimum and alternate in signs (starting with a negative) for H a maximum. 
The example problem in Appendix II shows how this condition can be satisfied 
for two control variables. 


The final necessary conditions to be discussed are concerned with the 
boundary conditions. The x.°, A. 0 , t 0 , and are chosen to satisfy the following 
relationships: 


F a (x°, t°) = 0 


YV 


i 

Z p 


dF 


a 


a-i 


a 


9x. 

1 


1 Id F \ n 

«<^Zp a M + z 

0=1 tn 1=1 


C. X. (to) = 0 


Gg(x f , t f ) - 0 


m 


\ <V 


= - c i - Z p 


’9G, 


W X i 


m 

H <V-Z p, 


9G, 




= 0 


The boundary conditions are actually only n + 1 independent relationships 
in x., A^, and t at t 0 and t^. because the multipliers p^ and p^ can be elimi- 
nated by using 1 of the n + l+i relations at t 0 to solve for the p^'s and m 

of the n + m + 1 relations at t- to solve for the p' s. When the p ’s and p 's 

f p a (3 

are eliminated, the necessary conditions to be satisfied by the optimization 
problem can be written in a more compact system of equations 
9H 


x i 9 A. 

l 


i = 1,2, • • • , n 


9H 

i 9 x. 

l 
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0 ** 


m 

s \ 


k = 1,2,. 


r 


n 

H a maximum for a minimization of S= Y c. x. (t ) 

. “ 11 f 

i=l 


and n 

H a minimum for a maximization of S = Yj C i X i ^ V 

i=l 

with boundary conditions 

Dj (x 0 , A 0 , Uq, t 0 ) = 0 j = i,2,--*,n+ 1 

Ej > Ujj t^) — 0 • 

A derivation of the preceeding necessary conditions is given in Appendix I. 


The next section of this report will be concerned with three of the numer- 
ical techniques that may be used to determine the x 0 , A 0 , t 0 , and t^ that will 
satisfy the 2n + 2 relations D. and E.. 


Methods for Solution 


Newton's Method with Numerical Derivatives, - The simplest and most 
straightforward method of satisfying the Dj and Ej relations is to evaluate 
numerically the effects of small changes in x.^, A.^, t 0 , and t^ on the and 

Ej relations. This is done by computing what is called a trial with estimated 

values for the unknowns x.^, A.„, t 0 , and t r . Naturally the D. and E. relations 

lO iO f j j 

will not be zero, but if the unknowns x.^ and A.„ are changed, small amounts 

iO iO 

one at a time, the changes can be used to determine the correct values for x. 

A. Q , t 0 , and t^. The initial trials and the 2n small changes in the unknowns 

allow a numerical determination of the partial derivative of D. and E. with re- 

i 3 

spect to x.„ and A.„. 

iO iO 

6 G • 6 • G 

To show how this can be done let D . , E. , D. , and E. at t n and t. be 

1 J . J .1 f 

the symbols used to denote the values of Dj, Ej, Dj, and Ej at t 0 and t^ in the 


iO’ 


Ax 


Ax 


trial computation. Then let lu and E ^ 6 * * * 10 , at t 0 and t^ be the symbols used 

to denote the values of D. and E at t n and t„ in the computation of the trial 

J J ° f 
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associated with the small change Ax 10 in the single unknown x 10 . This allows 
the following expressions to be written: 


9D . ( t 0 ) 


9x, 


10 


9E.(t f ) 


D AX 1° - D 6 

_J L 

Ax 10 

E AXl « - E 6 

_J 1_ 


9x 


10 


Ax 


10 


The other partial derivatives 


9D j (t 0 ) 


ox. 


iO 


9E . (t f ) 


9x 


10 


9D.(t 0 ) 
1 


9X 


iO 


9E . (t f ) 

and -rd® can 


0A. 


iO 


be determined by similar independent small changes in the x.^ and X.^. DTt 0 ) 
and E.(t 0 ) are also computed for the estimated values of the unknowns x.q. 


t 0 , and t f . 


Newton's iteration formula can then be used in the following form: 


0 = 


D 6 

1 


E 6 

J 


+ 


9D.(t 0 ) 

9E . ( t ) 
.1 { 

8x i0 


3 — — 


9X 


iO 


D® (t 0 ) 
0 




r 

0 


Ax 

6 ‘ ( v 

1 

$ 

At 0 

At f 




This formula is simply the linear terms of a Taylor Series expansion in the 
2n + 2 variables x.^, X.^, to, and t^. 

Solving the matrix equation for Ax.^, AX^q, A to, and At^ gives the 
corrections that are needed for another trial. 



AX 


iO 


At 

At 


9D.(t 0 ) 

J 


9x. 


iO 


9D . ( t 0 ) 


9X 


iO 


D®(t 0 ) 


9E . ( t ) 
3 f 

9X iO 


9E.(t { ) 


9X 


iO 


— 

-1 


0 


r D e l 

3 

. p 

E., V 
1 


E® 

3 J 
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The entire process can now be repeated with the preceding solution as new 
estimated values for x , X. q, t 0 , and t^. Convergence of the process to values 

of x.q, X.q, t 0 , and t^ that yields zero values for and can be greatly en- 
hanced by a simple modification. Instead of using zero in Newton's formula, 


the quantity K 


Atj become 


Ax i0 


D® 

J 

E 6 
L JJ 


is used. Then the expressions for Ax.^, AX i>n , At 0 , and 


iO’ 


A \o 

At 0 
At f J 


= - ( 1 - K) 


8D (to) 

A 


8D (t 0 ) 

J 


ax. 


10 


9E (t f ) 

J 


a x 


10 


9E (t f ) 

A 


D®(t 0 ) 


ax 


iO 


ax. 


10 


• p 

E. <V 


d: 


E 

L JJ 


The constant K must be chosen 

b e 

i 

to be zero unless the vector 


in the range 0 ^ K < 1. It is usually chosen to 
is not decreased on a particular iteration 


cycle. Then K can be increased as close to one as is necessary to assure that 


the computed values ofAx.„, AX.„, At. , and At. 

iO iO 0 f 


produce a decrease in 



Another modification to this procedure is called the Secant method. Its 
main advantage is that it does not require 2n integrations of the system of 
differential equations for each iteration cycle. After the first iteration cycle 
only one additional integration of the system of differential equations is used to 
modify the partial derivative matrix. A more complete explanation of this pro- 
cedure can be found in Reference 1. 


Newton's Method with Integrated Partial Derivatives . - The next method 
of solving the boundary value problem also uses Newton's iteration formula, 
but the partial derivative matrix is computed in a more accurate and efficient 
manner. A rigorous justification for the procedure is found in Reference 2, 
although most of the steps are intuitively obvious. The method involves chain 
rule differentiation of the boundary conditions. The system of equations which 
results is: 


10 


3D. (I®) 

nT 


9D.(t 0 ) 


9D.(t 0 ) 

I_ 


ax 


iO 


aE.(tj 

i_L 


n 


aD.(to) 


a^ ( t 0 ) 


)x i 


_ 9x io. 


3D. (to) 


ax. 


ox. 


iO 


8E.(tJ 

J_L 


9E. (t ) 
J l 


ax. 


iO 


ox l 


aE.(tj 
J f 


i° 

3Xj ( t f ) 


ax. 


ap.(to) 


ax. 


~ayt 0 f 

-f 

~0D (to) 


au R (t 0 ) 

8x i0 j 


L 9 \ . 


9x io_ 

ax, (t 0 f 


~aD.(t 0 )“ 


\<V 

_ 9 \o. 

■f 

9u k 


9 \o 


3x i 


_ * x i0 

axj(t f )' 


aE.(tj 

..J. f 


ax 


iO 




3E.(tJ 
J f 


ax 1 (t f ) 


d \ 


" x i0 


ax } ( t f )' 


3E.(t f ) 


ax 


iO 


9u k 


9E.(tJ 

] f 




1 = 1,2, • • • ,n 


The only unknowns that appear on the right-hand side of the preceding equations 

n r ~i r~ — » -» r~ — » 


are the matrices 
~d 

The matrices 


n 


9x i0 


ax. 


9x io 


\ 


9x i0 


and 


)U k 


9 \o 


,U k 


ax. 


iO 


,x i 


9 \o 


n 


ax. 


iO 


, and 


9 \ 


9 \o 


at t 0 and t^ 


can be determined at any time t 0 < t < t^. 
from chain rule differentiation of the third equation in the system of equations 
9H 


0 = 


)U k 


Differentiation yields: 



a 2 H 



[n] 

1 

a 2 H 


P»i 1 


a 2 H 


\»\] 


8u 9xJ 
m 1J 



+ 

^m 9 ^! 


no. 

+ 

9u 9u, 
m k 


_ 9x io_ 

r 

a 2 H 


r : 

5x l 1 


a 2 H 


rn i 


9 2 H 


M 

3u 9x, 

L m lj 


no. 

+ 

9U m 9 ^l 


no . 

+ 

9u 9u, 
m k 


n 0 . 


m = 1,2, • • r . 
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Therefore 



must still be determined. The matrices at t 0 are determined from the initial 
conditions. 


Y t ° ) =x io 

Y‘o> -\ 0 • 




Differentiation yields: 
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The matrices at t^ are determined by integrating the matrix differential 

equations obtained by differentiation of the first and second equations of the 
system of equations 


*1 

*1 


8H 

3X. 

1 

9H 

ax. 

1 


Differentiation yields: 


K1 

d 

9X 1 


a 2 H 

8x i0_ 

dt 

_ 8x i0. 


_ 8x i 8x i 

"ax’ ' 

d 

[ 8 N] 


a2H 

_ 8x io_ 

dt 

ax ioJ 


ax.axj 


iL 

r ax i" 


a2H. 

HsT 

0 

1 

dt 

! 

. 8 \oJ 


ax. ax, 

1 U 

a> 

j** 

i 

d_ 

~dx [ “ 


a 2 H 

ho. 

dt 

L ax ioJ 


ax.axj 


CP 1 
** 
j 


a 2 H 


K 1 

L ax i0. 

+ 

ax.ai^ 


_ 8x i0 


ax 


i°_ 

9x i 


a 2 H 


ax. 




,X 1 


ax 


iO 


a 2 H 


l8x i 8 “k 


8u k 


fex 


10 


a \o 


a z H 


ax. a 


i d \ 


au k 


9X io 


L* 

—1 

L. 


L 

-J 


CP 

h* 

1 1 


a 2 H 


" 8X 1 " 


a 2 H 

V 

» CP 
© 

L 


<r 

CO 


? X ‘o. 



Kj 


The initial conditions for the preceding system of matrix differential 


equations 


determined. Also, 


3^ (to) 


^(tof 


ax^to)" 

, and 


_ 8x i0 _ 

> 

dX 

L io J 

9 

_ 8x i0 

8X iO . 


' 8 \' 


8x i0 


and 


dX 


10 


have already been 


have already been determined for 


every to < t ^ t^; therefore the equations can be integrated to yield the desired 
results at t^. The desired results at t^ are values for the matrices 


ax 


10 


jV 

ax i0 


V 

ax i0 


and 


aX l 


ax io 


These matrices are substituted 
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9E.(t f ) 9E (t f ) 

into the formulas for — and . Then all of the elements are 

8x i0 8 No 

again available for the use of Newton's iteration formula in the following form 




9D.(t 0 ) 

9D.(t 0 ) 


r\ 


1 

D?" 

J 

~ Ax i0 _ 

A \o 

_ 

8x i0 

8 No 

U 


A t 0 
_ At f _ 


9E.(t.) 

y f 

»E <tj 

0 

• P 

E. <t f ) 




_ 8x io 

8 No 


Ei . 

J 


The computations for the partial derivatives can be performed along each 
trial trajectory, and the results are much more accurate than the numerical 
differentiation procedure used in the first method described for satisfying the 
boundary conditions. The modification to increase the chances for convergence 
suggested in the explanation of the first method is still very helpful with the new, 
more accurate partial derivatives. 


Modified Newton-Raphson Operator Method. -The final method to be 
discussed is a modification of the quasilinearization or Newton-Raphson operator 
method. The modified method is similar to the steepest descent procedure in 
that prespecified control functions are needed to start the iteration cycle, but 
the derivation is such that no backward integration is required. The first step 
in the explanation of the procedure is the linearization of the system of equations 
by using the first-order terms in a Taylor Series expansion in all the variables 
about the initial trial trajectory, 


• 

X. 

1 

7C. 

1 

0 

m 



9 2 H 

9X.9XJ 

(Ax 1 ) + 

9 2 H 

_ 8X i\ _ 

(A V 


9 2 H 

(Ax 1 ) - 

9 2 H 

(AAj) - 

9 2 H 

9x.9x 1 

ax.9^ 

ax.au^ 

9 2 H 

(AXj) + 

9 2 H 

(AAj) + 

9 2 H 

9u 9x, 
m 1 

9u 8A. 
m 1 

9u 9u, 
m K 


r 


(A V 

<A V 
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1 = 1, 2, • • • , n 

k = 1,2, •• •, r 


Ax i = *i - *? 

AXj = Xj - 

AU k = \ ' “k 

The last equation in the preceding expansions can be solved for Au^ and 

the results substituted into the first two equations. Then the first two equations 
can be arranged as a linear system of differential equations with time varying 
matrix coefficients: 

x. = [Aj (t)] x. + [B t (t)] A. + C u (t) 

= [ A 2 <t)] x. + [B 2 (t)] A. + C. 2 (t) . 

The coefficients A, B, and C are generated along the trial trajectory using the 
arbitrary guesses for x.^, A.^, to , t^, and the arbitrary control functions 

To obtain a general solution to the above linear system of differential equations 

P P 

one particular solution (x. and A. ) of the entire system must be generated 

by numerical integration and 2n particular solutions to the homogeneous part 

(i. e. , with the C., and C terms left off) must be generated. The initial 
ll 

conditions for the particular solution to the nonhomogeneous system are chosen 
to be the same as the initial conditions for the trial trajectory. The initial 
conditions for the 2n homogeneous solutions for convenience are chosen as 
follows and written in matrix notation: 

x (to)= [I] A <to)= [0] q = 1. 2, • • • ,n 

Also 

x. (V = [0] A (to) [I] s = 1, 2, • • • ,n 

IS lo 

Each of the columns in the previous set of matrices is considered to be a set of 
initial conditions to yield the 2n particular solutions to the homogeneous sys- 
tem. All of the numerical integrations can be performed at the same time, and 
the general solution is then written: 
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X.(t) 

+ 

II 

[V 0 ] 

x.(t) 

= x^ (t) + 

1 >> 1 
» M. 

l l 

The constants K 

and K 


K + 

q 


M 


K 


that the general solution will satisfy the and boundary conditions. To 

do this the boundary conditions are also expanded in a Taylor Series with only 
the first-order terms retained. This yields the following expressions: 


0 = D? 
J 


0 = E. 


3D. 


ax. 


r 9E® 1 



3D® 


3D 6 

Ax.(t 0 ) + 

- to 

ax. 

1 

AX. (t 0 ) + 
to 

3 

9 \_ 


Au k (t 0 ) + D®(t 0 )Ato 


'to 


ax. 

i 


Ax.(t f ) + 


r 8E® n 


ax. 

1 


AX^(t^) + 


K 


8e: 


9 \ 


Au k (t f ) + E® (t f )At f 


Ax.(t f ) = x.(t f ) - x®(t f ) 
AX.(t f ) = X.(t f ) -X®(t f ) 


Again the third equation in the system of equations is used to obtain 

terms of Ax. and AX.. Then the general solutions for x. (t) and 
ii i 

ated at t 0 and t^ are substituted into the Ax. and AX^ expressions, and the 

boundary conditions become 2n + 2 linear equations in the unknowns K , K , 
At 0 , andAtf, which are easily solved. 


au * ln 

X.(t) evalu- 
i 


Using the newly computed values for K , K , At 0 , and At f , the original 

C[ S X 

system of equation and the linear system with time varying coefficients are re- 
integrated. The initial conditions for the original system of equations remain 
the same, but the initial conditions for the linear system are given by evaluating 
the general solutions for x. (t) and X. (t) at t 0 . This yields 
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X.(t 0 ) 

i u new 

4-i 

n 

+ 

K 

q 

A. (to) 

i u new 

= *®(t<>) 


K 

s 

to 

u new 

= to 


At 0 

t, 

f new 

®JM 

II 

+ 

At f . 


During the reintegration process the third equation of the system of equations 
is again used to compute Au^, which is added to the old control function to 

-i _ . j.i. . a i jC a * r n. . *. a • 1 m.:- ~ 

pruuuut? tilt? ntsw uuntrui iuuuuuii iur me neAL a i*; Act nun uyu te* a hid piuucumc 

does not converge until- both the boundary conditions are satisfied and 

Oil 

= 0 over the interval t 0 < t £ t- 

8 \ r 

To increase the' chances of convergence the process can be made to creep 

8H 

toward a solution by not requiring that D., E., and — — be zero on a particular 

3 3 9 \ 

trial, but only that it be smaller than the previous trial (i. e., zero in the Taylor 


8H 

Series expansions for D., E., and - — is replaced by K 

3 3 dU k 

where 0 < K < 1 and 0 ^ k < 1) . This is the same idea that was suggested for 
the first and second methods of satisfying the boundary conditions. 


r D e i 

3 



and k 

'sh' 

E? 


9u k 

. 3. 




Computational Considerations 


All three of the preceding methods require that the units of length, mass, 
and time be adjusted or scaled such that each of the variables has the same order 
of magnitude. This scaling is needed mainly to retain good numerical precision 
when the matrix inverse in Newton’ s formula is taken or when the system of 
linear equations is solved in the third method. For trajectory optimization prob- 
lems this is easily accomplished as is shown in the Computational Procedure 
Section of Appendix IL 
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None of the methods is sensitive to the choice of the quantities that must 
be estimated to begin the iteration cycle, except the first method. If the Ax 
values estimated for the numerical determination of the partial derivatives are 
not chosen properly, this method may not converge. The second method has 
been found to converge for almost any values of the quantities that must be esti- 
mated, except all zeroes. For example a heliocentric transfer from the Earth’s 
orbit to the orbit of Mars converged in only 10 iterations with very crude guesses 
for the initial conditions. In Reference 3, this same problem with the same 
crude guesses took 13 iterations and the procedure used (which is similar to 
method three) is slightly more complicated to program on a computer. Because 
the second method has proven to be so effective, no effort has been made by the 
author to program the third procedure. The explanation of the third procedure 
is included as a generalization and extension of the ideas presented in Reference 
3, and also to point out this method's similarity to the steepest descent idea. 

Inequality constraints on the control variables can be handled very easily 
by all three methods. When the control variables are on the constraint boundary, 
H is not considered to be a function of the control U; therefore, all of the terms 
in the equations that contain first and higher order partial derivatives of H with 
respect to u are considered to be zero. To handle discontinuities or inequality 
constraints on the state variables or functions of the state variables and the 
control variables, modification of the necessary conditions is needed. The modi- 
fications have been developed in References 4, 5, and 6. All of the modifications 
to the necessary conditions usually require that additional boundary conditions 
be satisfied at other than the first and last points. Because the three methods 
discussed for satisfying boundary conditions have been formulated for two points 
(the first and last) , the extension to three or more points should be obvious. 

Also, because the iterations on the boundary conditions are performed at both 
ends of the trajectory, all three of the methods can be integrated backward or 
forward. The backward integration is helpful when the final conditions are 
extremely sensitive to changes in the initial conditions, but not "vice versa." 


CONCLUSIONS 


Practical methods for solving boundary value problems associated with 
the optimization of trajectories have been discussed. Actual experience with 
the construction of computer programs and the numerical results of computer 
programs has indicated that the second method described is usually the most 
effective for solving boundary value problems. Efficient use of the new IBM 
FORMAC computer language, which enables the computer to obtain functional 
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forms for the partial derivatives of functions with respect to the variables that 
appear in the function, will also assure that the second method becomes an 
easily programmed and economical computer program. More about the IBM 
FORMAC computer language can be found in Reference 7. 

The detailed application of the second method to a simple trajectory 
optimization problem is outlined in Appendix n. Because the theoretical formu- 
lation for more difficult problems is also available, subsequent efforts will be 
directed toward the application of the second method to these problems to gener- 
ate additional computer programs. 
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APPENDIX I 

DERIVATION OF THE NECESSARY CONDITIONS 
FOR THE OPTIMIZATION OF S 


The problem is to determine a u(t) that satisfies the system of differ- 
ential equations 


x. = f. (x, u,t) 
1 1 


with boundary conditions 


F a <*,,*,) = 0 

V x fV = 0 

and that maximizes or minimizes the quantity 
n 

s = E c i w • 

i=l 


i = 1,2,... ,n 

X = Xj.Xg,. • • ,x 

u = U^Ug,...,!! 


a = 1,2, • • • ,1 < n + 1 
/ 3 = i, 2 ,'*‘,m 2 =n+i 


Because the x. *s are assumed to be continuous, S can be rewritten in integral 
form. 1 


S 




c. x. 
x 1 


n 

dt + Z c n - x i 

i=l 


To examine the effects of the constraints on the maximization of S, a new 
function S' is defined. 



’ n n 1 

E c i x i + E \ - f i< x » u, t>] ; 

i=i i=i j 


dt 



A minimization of S' is equivalent to a minimization of S, if the boundary 
conditions and the system of differential equations are also satisfied. Assuming 
a solution to the problem statement exists, the optimum S' can be written as 
follows: 
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t f f n 


s ’ - 


Tj C. X. + y X. [ x. - f. (x,u,t)] 
f-' 11 l 1 i v 

i=l 1=1 


dt 


n _ i m 

+ Z c i + Z P a F a (xo»W + Z . 


i=l 


Q!=l 


0=1 


The bar over all of the variables means that they have their optimum value. 
Now a variation in all of the variables is written as follows: 

t. + At. f n 

r i f J v ,r . . v 
»•=»■+ ao' =i w. c. ix. + axj 

J- H. li l 

tg+Aio l i— i 


n 

+ Z (^- + AX.) [ (x. + Ax.) - f.(x + Ax, u + Au, t)] 

. j 1 1 1 11 


1=1 

n 


dt 


+ Z c i I MV + Ax (t 0 )] +Z (p +Ap ) F [x 0 (to)+Ax(t 0 ),t 0 +At 0 ] 

. 4 11 1 A QL 0 /. OL 


i=l 

m 


QL =1 


+ Z (Pp + A Pp> G^f^^V + Ax( V ’ *£ + A V 


In the preceding expression the variations for times greater than t^ and less 
than t 0 are taken from the values at t 0 and t^. For example 


Ax(t) = x(t) - x(t 0 ) 


t ^ t n 


Ax(t) = x(t) -x(t^) 


, etc. 


Subtracting S’ from S’ + AS' yields AS'. 


AS* 


t 0 f n n - 

f J Z c.(£ +Ax.) + Z (X. +AX.)[ (x. +Ax.) -f. (x + Ax, u 

J - < a * 1 . 4 i i i . LJ i i l i ii 

to+Ato I l— 1 1=1 


+Au,t)Udt 


f- 


n 


« ii 

Z c Ax + Z \ [Ax + f (x,u,t) - f (x + Ax, u + Au, t)] 

• i 1 A • i ^ * 1 1 


li=i 


i=l 


+ Z AX.[(x. + Ax.) - f.(x+ Ax,u + Au, t)]| dt 
i=i 
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t+At. n n _ _ I 

+ | i Y. c.(x. + Ax.) + V ( X.+AA ) [ ( x.+Ax.)-f. (x+Ax,u+Au,t )] >dt 

J 11 1 f-' 11 111 I 


i=l 


i=l 


n 


+ Yj C. Ax. (to) + 2 P y [x(t 0 )+Ax(t 0 ),t„ + At 0 ] - F [x(t 0 ) ,t 0 ] 
i=i 1 1 a = 1 “ 1 a 


+ Y AP a F a [x(t 0 ) + Ax(t 0 ),t 0 + At 0 ] 
a = 1 


m 


+ Y P B i G o[x(V + Ax(t f ), t 0 +At 0 ] - G [x (t f ), t f ] 
/3=i P ' P P 


m 

+ Y Ap rt iG n [x(t f ) + Ax(t # ), t # + At # ty 
/3=i 


"ft > 


If AS' is greater than zero for all the variations, S' will be a minimum and 
"vice versa. " To see how AS* can be made greater than or less than zero the 
following steps are taken. The expression for AS* is simplified by assuming 
that all of the variations considered are small. This allows Taylor Series 
expansions to be used, and all second-order terms can then be neglected. The 
following expansions are made with only linear terms of the Taylor Series re- 
tained. The symbol 5 means a smaller variation (meaning first-order terms 
are sufficient) than the symbol A . 


n 


f.(x + Sx, u + Au, t) = f. (x, u + Au, t) + Y 
1 1 1=1 


3F(x, u + Au,t) 


a 


*1 


n 8F a [x(t 0 ) t 0 ] 

F [ x ( t 0 ) + Sx( t 0 ) , t 0 + <5t 0 ] = F [ x ( t 0 ), t 0 ] + Y I ^ 

a a . =1 dx. 


a 


5x i 


<5x.(t 0 ) 



* 

- 

- 

/ 

_1_ 1 

n 

v 

aFjx(t 0 ),t 0 ] 

X.(t 0 ) 

fit 0 + 

+ 

L j 

. 1=1 

ax. 


i J 


- 


aF^fx (t 0 ) ,t 0 ] 


at 


n 


G [x(t f ) + 6x(t f ),t f t 6t f ] = G [x(t f ),t f ]+J < 


1=1 L 


f fito 


9G 


8x. 


6x.(t f ) 


r 

n 

J V 

3G [x(t f ),t f ] 

s<v 

r 6t^ + ' 

aG 0 (x(t f ),t f ] 

L 

i=l 

ax. 

at 

i 





fit. 


22 



Substitution into the expression for AS' and elimination of the second-order 
terms caused by the smaller variations (denoted by a 6) yield: 


n n 

6S' = -j J c x (to) + Yj \(t 0 ) ( x.(to) - f :.(xo,u 0 + Auo.toJlfsto 
1 i=l i=l 1 




t { f n 


XX XX p 

Z C. 6x. + Y X. 6x. + f.(x,u,t) - f.(x, u+AUjt) 
i=i i=l L 


n 8f. ( x, u+ Au, t) 

. y -i — 

u 8Xj 


l=i 

n 


<6x l j 


+ Yj 6 \- 1 K - f . (X, U + Au, t)] 
i=i 


r dt 


+ 1 


n n 

Z C iVV + Z VV [ x.(t f ) - f.<x f ,ii f + Au,t f )] l 6t f 
i=l i=l 


n 


+ Z c. 6x.(to) + Z P a * Z 

i=i a=l i=i 


(x 0 ,t 0 ) 


ax. 

i 


6x.(to) 


n 

i=l 

1 




dx. 


x.(t 0 )6t 0 + 


aF^(x 0 ,to) 


at 


<5tn > 


m 


+ Z^ 6 p a + Z^ 6 pp % f x rV 


a=i 
m 

+ Z PrA 


/?=! 


/3=1 


n 

Z 

i=l 


8 ViV 


aG B <vV 


at 


ax. 

1 


6t„ 


fix (tJ + Z 

i=l 


8G g( x f*V 


ax. 

1 


MV 6t f 


A combination of certain terms under the integral sign may be integrated by 
parts as is shown: 
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f 


" n 


" n 

Z (.C+XJ5X 

dt = 

V (C. + X.) 6x. 
Aj 1 l l 

i=l 


i=i 


f 

to 


t f 

n 

f 

y X. 6x. 
i i 

to 

1=1 


dt 


Substitution of this expression into 6S' and further rearrangement yields an 
expression of the form. 


6S' =< 


n__n__ _ __ 1 

Y c x (tj,) + 2 x.(t 0 )[x(t 0 ) - f (x 0 ,u 0 +Au 0 ,t 0 )] l 

i=l i=l J 

x. (t 0 )6t 0 + 


1 


♦ Z p„ Z 

a=i li=i 


9F a (x 0 ,t 0 ) 


ax. 

1 


aF a (x 0 ,to) 


at 


6t n 


6t« 


n 

z 

i=i 


\(to) + Y P 


l aF a (x 0 ,t 0 ) 


. a 9x. 

0=1 1 


Sx.(t 0 ) + Y 6p a F a (x..to) 

o=l 


3f (x,u + Au, t 


t rn n af (x,u + 

z 'VZ <v iT 

t 0 t i=1 L 1=1 1 

n 

+ Y ^.[t(x,u + Au,t) -f.(x,u,t)] 


6x. 

l 


i=l 

n 


- V 6\.[x. - f. (x, u + Au,t)] 
. lj a ii l 


dt 


i=l 


{ n n 

Z c i 5 i<V + z \ 

1=1 1=1 


(t f ) [x.(t f ) - f.(x f ,u f + Au f , t f ) ] s et f 


m ( n 

♦ z pjz 

(3=1 ' [ 1=1 


8G /3 (X f > V 


ax. 

1 


x.(t f ) 6t f + 


SG^Xf.tf) 


n 


z 

i=l 


- _ “ . dG B Cx rV 
°i + VV + 2 ^ to— 


0=1 




m 


St, 


<5x(tJ+ J 6 P« G « (x f‘V 

S=i p p 
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The only way that 6S’ can be made greater than or less than zero for either 
positive or negative values of the variations (the 6 values) is to make the 
coefficients of all the variations zero. This yields the following necessary 
conditions: 

F a (^t 0 > = 0 
_ _ 1 _ 9F (Xo,t 0 ) 

p„ ex. - 

a=l i 


n n 

V e v / T \ j. V _i_ a ». + n 

LJ • l_t ,v jV v l)' *jV^0>“0 ‘ •'O'J 

i=i i=l 




|y 

9F a (x 0 ,t 0 ) 

\h 

9x. 

l 

= 0 
n 

v fSei 


_ _ 9F (x 0 ,t 0 ) 1 

k i (t ° )+ ~^t~ 


> =0 


/3=1 P 


9X 

1 


n n 

Z o.k.(i f ) +Yj\ (9 [k i ( V " V^f’^f + Au f» V ] 

i=l i=l 


m 

♦ E p . 


If 


14 

9x. 

l 


VV * 


a( V x rV 

9t 


> = 0 


x. = f . ( x, u + Au, t) 

li ’ 


t n — t < t„ 


E \ 

1=1 


9ij ( x,u + Au, t) 


9x. 

l 


t fl < t < t. 


Then the expression for 6S’ becomes 


5S* f f |Z \ [f. (x, u + Au, t ) -f.(x, u, t))l 

t 0 li=l ] 


dt 
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Now SS’ will be positive if the variations Au make the variations in the 
differences of the f. either zero or negative. This means the variations in the 

differences of the f. must be negative for some finite interval in the time period 

t 0 :£ t :£ t^. This allows bounds to be placed on the control variables. From the 

opposite viewpoint, SS’ will be negative if the variations Au make the varia- 
tions in the differences of the 1 either zero or positive. The same considera- 
tions hold for bounds on the control variable. If AS f is positive for all small 
variations in u, then S' is a local minimum with respect to u and "vice 
versa. " Defining a new quantity H and combining and rearranging the previous 
results allows the final form of the necessary conditions to be written 

n 

H = E \ f i < X >M) 

i=l 

Then, 

x. 

1 

X. 

1 

Because H must be a maximum or a minimum with respect to u 
9H 

— = 0, and H a maximum for S a minimum 
8u 

H a minimum for S a maximum 


3H 

d\. 

l 

9H_ 

9x. * 

l 


F a ( x 0»t 0 ) = 0 


Y*°> “2,p a 


1 


a = 1 


9x. 


n 1 

E c n - *4 (to) + H(to> + E P 


9F ( x o»to) 


a 


i=i 

G^(x f ,t f ) = 0 


a=i 


a 9t 


= 0 


m 9G_(x.,tJ 


YV— -Sf- • h<«j) - p 


m 9G /3 ( x f , t f ) 


13 at 


= o 
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This development and form for the necessary conditions were chosen 
because they allow the part of the second variation associated with the control 
functions (i. e. , the maximization or minimization of H) to be examined and 
satisfied. In reality all of the second order terms in the Taylor Series ex- 
pansions should be retained and examined. Combinations of the second order 
terms would produce quadratic forms, and all of these would have to be positive 
or negative definite to assure a minimization or a maximization of S. From a 
computational standpoint, it is impossible to assure positive or negative definite- 
ness for any but the terms associated with the control functions alone. The 
other terms can be examined as a test, but there is no freedom to correct the 
terms if the test is not satisfied. Therefore, a solution is usually obtained 
satisfying the necessary conditions given, and physical reasoning is used to 
determine if the trajectory is acceptable without resorting to the extra effort 
necessary to test for the sufficiency conditions. 
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APPENDIX II 

APPLICATION OF THE NECESSARY CONDITIONS 
TO A TRAJECTORY PROBLEM 


A simple model for the equations of motion of a 
in three dimensions is given by the following system of 


.. F . MX 

x = — sin v - — 3- 
m *p R^ 


y = - cos Xp co S>Jr - fs 


zl = - — cos v sin v 

rv- 1 ''-n /v 'i 


m 


-JSL 

R 3 


rocket- powered vehicle 
differential equations: 


R = Vx? + y 2 + z 2 


m(t) = m 0 - m 0 (t-t 0 ) 
F ( t) = F o . 


These equations are based on the same assumptions that were used for the two- 
dimensional model described in the first part of the General Discussion. The 

control variables are x and x that locate the missile axis or thrust with 

P y 

respect to the coordinate system that is shown in the diagram below. 
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The following substitutions are made so that the system of differential 
equations will be in the form required for the application of the necessary con- 
ditions: 


Then 


Xj = x 

x 2 =y 

x 3 = z 
X 4 = X 

x s = y 

*€ = Z 

x 7 = m 
xe = F 
u l = X p 
*2=X y 

fl-A 

Xr 


SinUl " (x 4 2 +^+x 6 2 )% 


zr cos “■ cos Ui ■ 

t s = * cos u, sin a, - (Xt i l ^ + ^ > i/ t 

f 4 = Xj 

f 5 = x 2 
f 6 = x 3 

f 7 = -m 0 
f 8 = 0 


Now H can be written: 

H = * in +x *[{^) COB “i cos “» - 

+ [”(%) COS Ul Sin Uz ~ (X 4 Z 4 X 5 2 +X 6 2 ) 3 / 2 ] + ^ 4Xl + + ^6 X 3 - X 7 m 0 
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Then the IBM FORMAC language can be used to obtain fortran expressions for 
part of the necessary conditions as follows: 


8H 

X. = TT~ 

1 9A. 

1 

.. _ m_ 

i ~ 8x. 


i = 1 , 2, •• • ,8 


Expressions for the control variables that appear in the x^ and X^ equa- 
tions are obtained by solving the third equation in the system of equations for 
u* and u 2 , 

^2- = 0 = Aj p^ 8 ) cos u 2 - A 2 p^lsin u t cos U 2 + \ 3 f^ 8 -) sin U! sin u 2 
9uj \x 7 / \x 7 / \x 7 / 

|S_ _ o = - \ 2 cos Uj sin - A 3 cos u 2 cos u 2 

From the second equation 
A, 

tan iu = - -r* 

This is assuming that l"^ 8 ) and cos Uj are not zero, but L'Hospital's rule can be 

used to show that the expression is still true as both approach zero as a limit. 

From the previous relation 


sin u 2 


cos u 2 = - 


±\/ A2 2 + A 3 2 
+ A? 

’ W A2 2 + A3 2 


Substitution of these equations into the equation for — — gives: 

dUj 


At cos u t - f-— sin u t - sin u t = 0 

A2 A3 ±\j A2 + A3 

A t cos Uj ± N A 2 2 +A 3 2 sin u t = 0 


tan Uj = 


A2 2 + a 3 2 


*4 

SmUl ±n/ Xj 2 + Ag 2 + A 3 2 

iN / Ap 2 A 3 2 

COS U t = 

±4 Aj 2 + A2 2 + A3 2 


The ambiguities on the signs of sin and cos of u t and u 2 can be resolved by 


9 2 H 


9 2 H 


Because it is desired to minimize 


examining — 2 and — y 
8 

3 = c t x ,- v V where e t = c 2 - e 3 = e 4 - c 5 - e 6 - e 8 - 0 
i=l 


< 3 2 H __ j 

dLl 1U <J* — — ± 9 ~ 9 ttiiU 


’ 9UJ 2 


a 2 H 


2 must be negative. The above statements say that S = - m(t.) is to be a 

auz 1 

minimum. This is equivalent to stating that m(t^) is to be a maximum. For 

- m(t,) to be a minimum H during the period t^ ^ t < t f must be a maximum; 

.. , a 2 H , a 2 H 

therefore - — and — — 

3 m au 2 

tn — t — t 


must be negative at all times in the interval 


f ' 

3 2 H 

au, 2 


= - A 2 |^*| cos Uj cos u 2 + A 3 cos Uj sin u 2 


x 7 


= (^ L ) ^ Sin U 1 + C0S U 1 [ “ *2 COS u 2 + *3 sin u 20 


= “ *2 (“jjj) COS U 1 COS U 2 + *3 ( j^) COS u i sin 


u 2 


- 2k 


{ cos Uj [ - A^ cos u 2 + A 3 sin u 2 ]} 


From the physics of the problem it is obvious that will always be positive. 

Also, the positive sign is arbitrarily chosen for cos u t . 
n/A , 2 + A3 2 

COS Uj = 


4 A! 2 + a 2 2 + a 3 2 


a 2 H 


Then — 2 will be negative if the following signs are chosen for sin u 2 and 

3u 2 

cos u 2 
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cos u 2 = 


^2__ 

^ x 2 2 + x 3 2 


sin u 2 = 
8 2 H 


4 X 2 2 + X 3 2 

Also, -^ 7 y will be negative if the following signs are chosen for sin u 4 . 

h 


auj' 

sin Uj = 


4 Xj 2 + \ 2 + X 3 2 


9 2 H 


a 2 H 


The preceeding choice of signs make Qu 9 and 9 9u zero; 

1 u 2 2 1 


9 2 H 


9 2 H 


therefore H is a maximum because ~ — T~ and — j— are both negative. With 

8Ul 9 u 2 

the expressions for the control angles defined, the boundary conditions can now 
be discussed. 


For this problem the boundary conditions at t 0 are assumed to have the 
form shown below: 


F i = x 4 - X(, = 0 

f 2 = * 2 - y 0 = 0 
F 3 = x 3 - z 0 = 0 
F 4 = x 4 - Xq = 0 
f 5 = x 5 - y 0 = 0 
F 6 = Xg - z 0 = ° 

F 7 = x 7 - m 0 = 0 
f 8 = x 8 - F 0 = 0 
F 9 = t - t 0 = 0 

Then the transversality conditions become 


Xi(t 0 ) = P! 

X 2 ( tn) = po 


Xsl t 0 ) - p 5 

^■e(W = Ps 
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^3(^0) ~ P3 
A. 4 (t 0 )= P 4 

m 0 + H(t 0 ) + Pg = o 


fy(to) ~ 1*7 
^(t-o) = Rb 


The last equation determines pg , and the other eight p's or Vs must be 
determined to satisfy the boundary conditions and transversal! ty conditions at 
tj. For this problem the boundary conditions are assumed to have the form 

Gj = X4 2 + X5 2 + Xg 2 - R 2 = 0 

f 

rj_ = y.2 4- v~2 + v„2 — v 2 _ n 

~ l *-1 • --a * f 

f 

G 3 = X A X4 + X 2 X5 + x 3 Xg - (Rv cos',?) = 0 

t f 

The quantities R 2 , v 2 , and (Rv cos 1?) are desired constants that charac- 
f • f t 

terize a certain orbit where R is the radius, v is the velocity, and Rv cos ■& 

is the dot product of R and v . Then the transversality conditions become 

Mt f ) = 2P2*i + P3X4 

A 2 (t f ) = 2P2X2 + p 3 x 5 

X 3 (t f ) = 2 p 2 x 3 + p 3 xg 

^ 4 (t f )= 2p i x i + p 3 x 1 

X 5 (t f ) = 2 Pl X5 + p 3 x 2 

\(t f ) = 2 Pl Xg + p 3 x 3 

Mt f )= 1 

A 8 (t f )= 0 
H (t f )= 0 

This is nine transversality conditions and three boundary conditions for 
a total of twelve conditions that must be satisfied by the eight initial Vs or p’s 
and the three p*s at t^ and t f . The three p’s at t^ can be eliminated by 

using three of the transversality conditions to solve for the p's in terms of 
x's and X's. To do this in an easy manner the following steps are taken: 
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*2 = 2 P2 *2 + P3 x 5 


A-5 = 2 P 1 X5 + P3 x 2 


Vector algebra is used on these two vector equations to eliminate the 

three p's. Crossing the first vector equation with x 2 yields 

x 3 


x 2 X *2 = 2p 2 X 2 X x 2 + P 3 x 2 X x 5 


= P3 *2 x x 5 


Crossing the second vector equation with xf yields 


2pi xj x x 5 + p 3 xg x xj. 


= P3 l x 5 1 x | x 2 


Adding the two resulting equations gives 
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X 1 


V 


*4 


V 

x 2 

X 

*2 

+ 

Xg 

X 

*5 

x 3 


^3 




A 



V 


x 4 _ 




V 

= P3 

x 2 

X 

Xg 

+ P 3 

x 5 

X 



x 3 


X€ 


xe_ 


x 3_ 


Performing the cross multiplications yields the three desired independent 
scalar equations: 


X 3 X 2 - \ 2 X 3 + XgXg - XgXg = 0 
- X 3 x t + Xjx 3 - Xgx 4 + X^Xg = 0 
X 2 x i - X^ + X5X4 - X 4 xg = 0 


Now the eight X's or p's at t 0 and t^ are used to satisfy the follow- 
ing nine conditions at t^. 

E 4 = x 4 2 + x§ 2 + Xg 2 - R 2 = 0 

l f 

E 2 = x 4 2 + Xj 2 + x 3 2 - v 2 =0 

f 

E 3 = x 4 x^ + x 2 Xg + x 3 Xg - (Rv cos i? ) = 0 

t f 

E 4 = X 3 x 2 - X 2X 3 + Xgx 5 - XgXg = 0 

E 5 = - X 3 X| + X t x 3 - AgX 4 + X^tXg = 0 

Eg = ^2 x i “ ^1*2 + ^-5 x 4 ~ ^4 X 5 = ® 

E 7 — Xy — 1 — 0 


Eg - Xg 
Eg = H = Xj 


t) 8i 


+ A.O 


ft) 


sin Uj - 


COS Uj cos u 2 - 


M Xa 

(x 4 2 +x 5 2 +x € 2 ) 3 / 2 
M Xr 

(X 4 2 +X 5 2 +Xg 2 ) 3 / 2 
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+ A. 4 X 4 + A. 5 X 2 +A 6 x 3 -A 7 m 0 = 0. 


+ \ 3 



cos u 4 sin u 2 


M gfi ' 

(X 4 2 +X 5 2 +X6 2 ) 3 / 2 


The preceeding discussion and the discussion in the main part of this 
report should be sufficient for an understanding of the following computational 
procedure that was used in constructing the actual computer program. 

Computational Procedure 

Preload 


x^ = x (m/sec) 

V 


x 2 ° = y (m/sec) 

*2° 


x 3 ° = z (m/sec) 

*3° 


x 4 ° = x (m) 

\ 4 ° 


C71 

O 

II 

? 

*5° 


Xg 0 = z (m) 

*6° 


n , kg sec 2 x 

x 7 u - m ( ) 

1 m 

\ 7 ° 


xg 0 = F (kg) 

>* 

CO 

o 


F o = (kg/ sec) 

V 

cutoff altitude squared (m 2 ) 

. , kg sec , 

m o - — ) 

m 

f 


V 

f 

cutoff velocity squared (m 2 /sec 2 ) 

GM (m 3 /sec 2 ) 

(Rv 

cos t?) (determines path angle at 



f cutoff) 

t 0 (sec) 

K 



t f (sec) 

Tolerance = . 5 x 10~ 6 


At (sec) for integration 
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Preload Computations 


All input must be scaled 



Kj = 1.5698587 X 10~ 7 

(scales length) 

- 

K 2 = 1.241825 X 10" 3 

(scales time) 


K 3 = 0. 16001332 x 10 -4 

( scales weight) 

1 

1 

The proceeding scale factors are for near earth trajectories, and they 
cause the initial radius ( \l x 4 2 +x 5 2 +x € 2 ) , the initial mass, and |i to be unity for 
the test case data given at the end of this computational procedure. 



V = Xg 0 (Kg) 
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GM 

= GM 

(£) 

to 

= to 

(K 2 ) 

tf 

= *f 

(K 2 ) 

At 

= At 

(K 2 ) 

R t 2 

f 

-s 

(K! 2 ) 

V t 2 

f 

■v 

m 

(Rv 

COS i?) 

*. ■ (i 



Preload Computations for the Isolation Routine 


Set up the following matrices: 


9x. 


ax. 


ax. 

i_ 

ax.° 

i 


= [0] 


i (number of rows) = 8 
j (number of columns) = 8 


= [I] 


Suggested Order for Computing Line "n" for the Isolation Routine 


Y\ ■ 


~8x.~ 


ax. 

l 

, and 

l 

ax.° 

ax.° 

. j J 


L j _ 


are brought forward from the previous (t) 


to a new ft + At) by Runge-Kutta integration. 
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The following equation is to be used with the IBM FORMAC language to 
obtain functional (fortran) expressions for the partial derivatives needed for 
the subsequent calculations: 


H = X 


1 (^) Sin Ul ' (x 


+ ^ 


(a-) c°s u j cos 


+ a, -jj*- cos u, sin u 2 - ( Aa/z + 


(1) sm u » " ,/ , ' i - , ii . , i 

v ''Z 3 

\/ X,^ + X 3 2 


(2) cos u 4 = 

(3) sin U 2 = 


si X t 2 + X2 2 + X3 2 

— -_*a 

n/Xz 2 + X 3 2 


2^ 


(4) cos u, - 


... • 8H 

< 5) x i = W. 

1 


(6) x; 


9H_ 

3x. 

1 


i = 1,2, •**,8 
i = 1,2, •••,8 


Construct the following matrices with the IBM FORMAC language: 


3 2 H I 

8 \ 8x iJ 

a^H 1 

5 V\ J 


k(rows) = 1,2 


l(columns) =* 1,2, • • • , 8 


8u_ 8u 
k mj 


k(rows) = 1,2, 
m(columns) =1,2 
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i(rows) = 1,2, • • • , 8 
j ( columns) = 1 , 2 , • • • , 8 



tion is stopped when t =* . 


All of the following equations are to be computed at the time t - t^ on a 
particular trial: 

(1) Ej = X 4 1 2 3 4 5 + x 5 2 + x 6 2 - R 2 

f 

(2) E 2 = Xl 2 + x 2 2 + x 3 2 - v 2 

f 

(3) E 3 = x 4 x 4 + X 2 X 5 + x 3 x 6 - (Rv cos 

(4) E 4 = \ 3 x 2 - AgXjj + XgX 5 - A 5 Xg 

(5) E 5 = - \ 3 Xi + Ajx 3 - X 6 x 4 + \ 4 Xg 
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(6) Eg = A 2 x 1 - X t x 2 + \ 5 X 4 - \ 4 x 5 

(7) E 7 = A 7 - 1 

(8) E 8 = 

< 9 > E S = *,[(*) 

+ ** [(^) cos “i cos “« - 

+ x > [if ,) cos “• sln '“* ■ (4ww> l/2 

m 

+ A 4 x 7 + X 5 x 2 + AgXg - A 7 m 0 + X 8 F 0 
(10) Define E = [Ej, E 2 , E 3 , E 4 , E 5 , Eg, E 7 , Eg, Eg] 

Construct the following matrices or vectors using the automatic partial 
differentiation routine: 



R(rows) = 1, 2, * * • , 9 
i (columns) = 1,2,* #, ,8 
j (columns) = i,2,***,8 
k (columns) = 1,2 

Then compute: 



( 12) [w j = f(A 2 sin u t sin u 2 + \ 3 sin Uj cos u 2 ) 

I ^ m J 1 ( \ 3 cos u 7 sin u 2 - sin u 7 - Ag cos u 7 cos u 2 ) 

(A 3 cos Uj sin u 2 - A 2 cos u 7 cos u 2 ) 
(Ag sin u t sin u 2 + A 3 sin Uj cos u 2 ) 
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(13) ^ = 


w 


km 


-l 


X 2 cos Uj sin u 2 + X 3 cos Uj cos u 2 

X 2 sin u 4 cos u 2 - \\ cos Uj - X 3 sin Uj sin u 2 

3E, 



K 1 

• 

K 1 

• 

[«ki 

11 

w 

9x. 

1 

X. + 
1 

_ 9 V 

\ + 

9X. 

L 1 J 


X. + 


R 


i 9 t 


(15) A E = E 

Compute | AE| = \Jfef + E 2 ^ + E 3 2 + E 4 2 + E 5 2 + E e 2 + E 7 2 + Eg 2 + E 9 2 

(16) 


<7 


9Et} 

-K P 

Zj 

qs 


9 X.® ’ R 

L J 


L j J 


If |AE) is < tolerance go to converged case run saving X.°(old) and 
tj. (old). If not continue: 1 


(17) P = 


qs 


-1 — 
E 


Compute z^J x in double precision 


-1 . 


[ Z qs] [ Z qsJ 


-1 


and test fz 1 fz 1 =[!]. 


where 


p — [AXj®, AX 2 ®, AX 3 °, AX 4 ®, AX 5 ®, AXg®, AX 2 ®, AXg®, At^] . 


On the first trial skip the next three "if" statements and go to "now for 
a new trial" : 


I AE 


If 


n 


I AE 


n -1 


l A E 

If — = 


n 


I AE 


> (1-K+y K) set K = y~ 8 . 


< ( 1-K+ — K) set K = 1. 8 K . 
b 


n -1 

If K > 1 set K = 1 


Now for a new trial: 

X^° (new) = Xj° (old) - K A X. 
At f (new) = t f (old) -KAt f 
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Converged Case Run 


Reintegrate the trajectory with the converged A. (t^ and t and compute 
additionally at each At . 1 

R = n/x 4 2 + x 5 2 + xg 2 
v = n/x , 2 + x 2 2 + x 3 2 


a _ x i x< + x ? x 5 + x^ 
Rv 


sin t? = \ll - cos 2 1 ? 

■ , / sin $ \ 

i cos £) 

, 1 sin u< \ 

u< = arc tan 1 L ; 

V cos Uj / 

, ; sin u, \ 

u, = arc tan { L 1 

. cos u 2 / 

Print out at each At: 


(convert to degrees) 
(convert to degrees) 
(convert to degrees) 


(1) On the first step of each trial print 


t, x. , x.,A.,X. , 
111 1 


9u 


[3x. 1 


[ax. ' 


‘ax. ' 


raxM 

k 


1 


1 


1 

, and 

1 

ay 


0A.° 

J 

ax.° 

> 

ax.° 

ax.° 


J _ 


i 


j _ 


j j 


(2) At t = tj. on each trial print all of (1) plus the following 



1 

Qt 

w 

CO 

1 

• 

r -1 


r -j 

-1 

r -j 


r n 

E , AE , 

ax.° 
j _ 

’ V e r' 

Z 

qs 

j 

Z 

qs 

> 

Z 

qs 

* 


Z 

qs 


, and 


(3) At each At of the converged case run print 


t, u*, u 2 , R, v, J, x., A.,x., X. 
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1 a 


Data for an Example Trajectory 


to - 

150. 01366 


X 1 = 

2378.9375 


x 2 = 

1181.8890 


x 3 = 

-1259. 8308 


X4 = 

143391. 19 


x 5 = 

6447134. 8 


xe = 

-40967. 611 


x 7 = 

6361. 8535 


x 8 = 

40599. 685 


F 0 = 

0.0 


m 0 = 

9. 6776706 


GM- 

. 39860160 

x 10 15 

V = 

1.0 


II 

O 

4 * 

1.0 


1! 

CO 

<< 

1.0 


X4 = 

1.0 


^5 = 

1.0 


^6 = 

1.0 


x 7 = 

1.0 


*8 = 

1.0 


11 

Jr 

650 


At = 

1.0 


Tolerance = . 5 x 10 6 


R t 2 

= .44193245 x 10 14 

1 



V t 2 

= , 59878362 

X 

O 

OO 

1 



(Rv cos 1 ?) = 0 





K = . 1 
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